###################################################################
###Replication: The Politics of International Peace and Security###
###################################################################

#Load and/or install required packages#
if(!require(tidyverse)) install.packages("tidyverse")
library(tidyverse)

if(!require(MASS)) install.packages("MASS")
library(MASS)

###--load data at the UNSC-SB level of analysis. --###
#NOTE: Instructions to load data at the UNSC meeting level below.
#NOTE: save data to file and set save location as working directory
sc_subs <- read.csv("unsc_sb-level_final.csv")

###--Transform to use only UNSC-SB level data without yearly aggregate data--###
#sc_subs <- sc_subs %>% dplyr::select(Case:Notes)

###--Load UNSC meeting level data--###
#NOTE:This level is used for one robustness check in ISQ appendix (table 13)
#meeting_level <- read.csv(file = "unsc_sb-meetinglevel_final.csv")

###--Start ISQ Replication--##
###--For ISQ replication use UNSC-SB aggregate level data--###
sc_subs_yearly <- sc_subs %>% dplyr::select(YearProposed, sub_count_total:low_strength_created) %>%
                  group_by(YearProposed) %>%
                  filter(row_number()==1)


###--Figure 1--###

#transform data for graphing
sub_proportion_yearly <- sc_subs_yearly %>% dplyr::select(c(YearProposed, sb_proposed_sum, sb_created_sum)) %>% 
  pivot_longer(-YearProposed, names_to = "count_subs", values_to = "count")

#code for figure 1
ggplot(data=subset(sub_proportion_yearly), aes(x=YearProposed, y=count, fill=count_subs)) + geom_bar(stat="identity", position = position_stack(reverse=TRUE)) +
  scale_x_continuous(breaks = seq(1970, 2020, by=5)) + scale_y_continuous(breaks = seq(0, 14, by=2)) +
  scale_fill_manual(labels = c("Created", "Total Proposed"), values = c("darkblue", "steelblue")) + theme(legend.position = "bottom", axis.text.x=element_text(angle=90)) + 
  guides(fill=guide_legend(title=NULL)) + ylab("Subsidiary Bodies per Year") + xlab("Year") + theme_bw()


###--Figure 2--###

#transform data for graphing
sub_proportion_yearly_str <- sc_subs_yearly %>% dplyr::select(c(YearProposed, high_strength_created, med_strength_created, low_strength_created)) %>% 
  pivot_longer(-YearProposed, names_to = "count_subs", values_to = "count")


ggplot(data=subset(sub_proportion_yearly_str), aes(x=YearProposed, y=count, group=count_subs)) + geom_line(aes(color=count_subs), size = 1.15) +
  scale_x_continuous(breaks = seq(1970, 2020, by=5)) + scale_y_continuous(breaks = seq(0, 14, by=2)) +
  scale_color_manual(breaks = c("high_strength_created", "med_strength_created", "low_strength_created"), 
                     labels = c("High", "Med", "Low"), values = c("darkblue", "steelblue", "skyblue"),
                     guide_legend(title = "")) +
  ylab("Subsidiary Bodies per Year") + xlab("Year") + theme_bw()


##--Main Analysis--##

##--Table 4: SB proposal--##
model1 <- glm(sub_count_total ~ p5_avg + meeting_number + veto_sum + YearProposed, family="poisson", data = sc_subs_yearly)
summary(model1)
logLik(model1)
AIC(model1)

model2 <- glm(sub_count_total ~ p5_avg + meeting_number + veto_sum + YearProposed + log(pk_forces) + cold_war_dum, family="poisson", data = sc_subs_yearly)
summary(model2)
logLik(model2)
AIC(model2)

model3 <- glm(sub_count_total ~ us_rus_dist + meeting_number + veto_sum + YearProposed, family="poisson", data = sc_subs_yearly)
summary(model3)
logLik(model3)
AIC(model3)

model4 <- glm(sub_count_total ~ us_rus_dist + meeting_number + veto_sum  + YearProposed + log(pk_forces) + cold_war_dum, family="poisson", data = sc_subs_yearly)
summary(model4)
logLik(model4)
AIC(model4)

##--Table 5: SB Creation--##
model5 <- glm(sb_created_sum ~ p5_avg + meeting_number + veto_sum + YearProposed, family="poisson", data = sc_subs_yearly)
summary(model5)
logLik(model5)
AIC(model5)

model6 <- glm(sb_created_sum ~ p5_avg + meeting_number + veto_sum + YearProposed + log(pk_forces) + cold_war_dum, family="poisson", data = sc_subs_yearly)
summary(model6)
logLik(model6)
AIC(model6)

model7 <- glm(sb_created_sum ~ us_rus_dist + meeting_number + veto_sum + YearProposed, family="poisson", data = sc_subs_yearly)
summary(model7)
logLik(model7)
AIC(model7)

model8 <- glm(sb_created_sum ~ us_rus_dist + meeting_number + veto_sum  + YearProposed + log(pk_forces) + cold_war_dum, family="poisson", data = sc_subs_yearly)
summary(model8)
logLik(model8)
AIC(model8)

##--Figure 3--##
##use model 3 for expected SBs proposed##

#create prediction frame
max_dist <- max(sc_subs_yearly$us_rus_dist)
min_dist <- min(sc_subs_yearly$us_rus_dist)

s1 <- data.frame(us_rus_dist = seq(from = min_dist, to = max_dist, by = 0.25),
                 meeting_number = mean(sc_subs_yearly$meeting_number),
                 veto_sum = mean (sc_subs_yearly$veto_sum),
                 YearProposed = mean(sc_subs_yearly$YearProposed))

predictions <- predict(model3, s1, type="response", se.fit=TRUE)

pred_df <- data.frame(us_rus_dist = seq(from = min_dist, to = max_dist, by = 0.25),
                      count = predictions$fit, 
                      lower = predictions$fit - (1.96*predictions$se.fit),
                      upper = predictions$fit + (1.96*predictions$se.fit))

#graph the predicted counts
ggplot(data = pred_df, aes(y=count, x=us_rus_dist)) + 
  geom_line() +
  geom_ribbon(aes(ymin= lower, ymax= upper), fill = "midnightblue", alpha = 0.25) +
  scale_x_continuous(breaks = c(2,2.5,3,3.5,4,4.5)) +
  scale_y_continuous(breaks= c(2,4,6,8,10,12)) +
  xlab("US-Russia Preference Distance") +
  ylab("Expected Subsidiary Bodies") + theme_bw()

##--Table 7: SB proposed by body strength --##

models1 <- glm(high_strength_created ~ us_rus_dist + meeting_number + veto_sum + YearProposed, family="poisson", data = sc_subs_yearly)
summary(models1)

models2 <- glm(med_strength_created ~ us_rus_dist + meeting_number + veto_sum + YearProposed, family="poisson", data = sc_subs_yearly)
summary(models2)

models3 <- glm(low_strength_created ~ us_rus_dist + meeting_number + veto_sum + YearProposed, family="poisson", data = sc_subs_yearly)
summary(models3)


##--Supplemental Appendix--##

##--Table 1 appendix: SB level descriptive stats--#
sc_subs_descriptives <- sc_subs %>% dplyr::select(YearProposed, YearEstablished, YearTerminated, Resolution,
                                                  Nay.Votes, Abstentions, p.5.Vote, US.Veto, Russia.Veto, Strength)

summary(sc_subs_descriptives)
sapply(sc_subs_descriptives, sd, na.rm=T)

##--Table 2: SB type 2--##
summary(as.factor(sc_subs$Subsidiary.bodytype2))

##--Table 3: SB type 3--##
summary(as.factor(sc_subs$Subsidiary.bodytype3))

##--Table 4: SB purpose--##
summary(as.factor(sc_subs$Purpose))

##--Table 5: SB proposed and created by decade--#
creation_rate_decades <- sc_subs %>% mutate(decade = dplyr::case_when(
                                      YearProposed <= 1979 ~ '1972-79',
                                      YearProposed > 1979 & YearProposed <= 1989 ~ '1980-89',
                                      YearProposed > 1989 & YearProposed <= 1999 ~ '1990-99',
                                      YearProposed > 1999 & YearProposed <= 2009 ~ '2000-09',
                                      YearProposed > 2009 & YearProposed <= 2020 ~ '2010-20'),
                                      decade = as.factor(decade)) %>%
                                      group_by(decade) %>%
                                      summarize(proposed = n(),
                                            created = sum(sb_created_dum),
                                            success_rate = created/proposed)

##--table 6: Meetings and regional focus--##
#NOTE: Bring in UNSC Meeting level data

meeting_level <- read.csv(file = "unsc_sb-meetinglevel_final.csv")
summary(as.factor(meeting_level$Region))

##---Figure 2 appendix: Distribution of Proposed Subsidiary Body by Country of Interest
if(!require(maps)) install.packages("maps")
library(maps)

#format data for the plot
df1<- data.frame(table(sc_subs$Country.1))
df2<-data.frame(table(sc_subs$Country.2))
df3<-merge(df1,df2, by ="Var1", all=T)
df3$Freq.y <-replace(df3$Freq.y, is.na(df3$Freq.y), 0)
df3$Freq.x <-replace(df3$Freq.x, is.na(df3$Freq.x), 0)
df3$count<-df3$Freq.x+df3$Freq.y
colnames(df3)[1]<- "Country"
df3$Country <- recode(df3$Country, "United Kingdom" = "UK") #United Kingdom and United States are not recognized by the program, recoded to acronyms 
df3$Country <- recode(df3$Country, "United States" = "USA")
df3 <- df3[-1, ]
world_map <- map_data("world")
world_map <- subset(world_map, region != "Antarctica")

#create the plot
ggplot(df3) +
  geom_map(
    dat = world_map, map = world_map, aes(map_id = region),
    fill = "white", color = "#7f7f7f", size = 0.25
  ) +
  geom_map(map = world_map, aes(map_id = Country, fill = count), size = 0.25) +
  scale_fill_gradient(low = "#efc9cb", high = "#a20203", name = "Count") +
  expand_limits(x = world_map$long, y = world_map$lat)+theme_void()

##--Table 7 appendix: bodies proposed by region--##
proposal_region <- sc_subs %>% mutate(Region = as.factor(Region)) %>%
                    group_by(Region) %>%
                    summarize(proposed_count = n())

#--table 8 appendix: descriptive stats
descriptives <- subset(sc_subs_yearly, select = c(sub_count_total, sb_created_sum, sb_proposed_sum, p5_avg, us_rus_dist, meeting_number,
                                                  veto_sum, cold_war_dum, pk_forces, high_strength_proposed, med_strength_proposed, low_strength_created))
summary(descriptives)
sapply(descriptives, sd, na.rm=T)

#stargazer(as.data.frame(descriptives), type = "html", digits=2, out = "descriptives_manuscript_analysis.html",
#          covariate.labels = c("SBs Total", "SBs Created", "SBs Not Created", "P-5 Distance", 
#                               "US Russia Distances", "SC Meetings", "Vetoes", "Cold War", "PK Forces", "High Strength", "Medium Strength",
#                               "Low Strength"))

##--Table 9 appendix: SBs proposed using nbreg -- NBREG
library(MASS)

modelnb1 <- glm.nb(sub_count_total ~ us_rus_dist + meeting_number + veto_sum + YearProposed, data = sc_subs_yearly)
summary(modelnb1)
logLik(modelnb1)
AIC(modelnb1)

modelnb2 <- glm.nb(sub_count_total ~ us_rus_dist + meeting_number + veto_sum  + YearProposed + log(pk_forces) + cold_war_dum, data = sc_subs_yearly)
summary(modelnb2)
logLik(modelnb2)
AIC(modelnb2)

modelnb3 <- glm.nb(sub_count_total ~ p5_avg + meeting_number + veto_sum + YearProposed, data = sc_subs_yearly)
summary(modelnb3)
logLik(modelnb3)
AIC(modelnb3)

modelnb4 <- glm.nb(sub_count_total ~ p5_avg + meeting_number + veto_sum  + YearProposed + log(pk_forces) + cold_war_dum, data = sc_subs_yearly)
summary(modelnb4)
logLik(modelnb4)
AIC(modelnb4)

##--Table 10: SBs Created using NBREG--##
modelnb5 <- glm.nb(sb_created_sum ~ us_rus_dist + meeting_number + veto_sum + YearProposed, data = sc_subs_yearly)
summary(modelnb5)
logLik(modelnb5)
AIC(modelnb5)

modelnb6 <- glm.nb(sb_created_sum ~ us_rus_dist + meeting_number + veto_sum  + YearProposed + log(pk_forces) + cold_war_dum, data = sc_subs_yearly)
summary(modelnb6)
logLik(modelnb6)
AIC(modelnb6)

modelnb7 <- glm.nb(sb_created_sum ~ p5_avg + meeting_number + veto_sum + YearProposed, data = sc_subs_yearly)
summary(modelnb7)
logLik(modelnb7)
AIC(modelnb7)

modelnb8 <- glm.nb(sb_created_sum ~ p5_avg + meeting_number + veto_sum  + YearProposed + log(pk_forces) + cold_war_dum, data = sc_subs_yearly)
summary(modelnb8)
logLik(modelnb8)
AIC(modelnb8)

##--Table 11 appendix: SB Creation, OLS--##
model9 <- lm(sb_created_sum ~ us_rus_dist + meeting_number + veto_sum + YearProposed, data = sc_subs_yearly)
summary(model9)

model10 <- lm(sb_created_sum ~ us_rus_dist + meeting_number + veto_sum  + YearProposed + log(pk_forces) + cold_war_dum, data = sc_subs_yearly)
summary(model10)

model11 <- lm(sb_created_sum ~ p5_avg + meeting_number + veto_sum + YearProposed, data = sc_subs_yearly)
summary(model11)

model12 <- lm(sb_created_sum ~ p5_avg + meeting_number + veto_sum  + YearProposed + log(pk_forces) + cold_war_dum, data = sc_subs_yearly)
summary(model12)

##--Table 12 appendix: SB strength, proposed

model13 <- glm(high_strength_proposed ~ us_rus_dist + meeting_number + veto_sum + YearProposed, family = "poisson", data = sc_subs_yearly)
summary(model13)
logLik(model13)
AIC(model13)

model14 <- glm(med_strength_proposed ~ us_rus_dist + meeting_number + veto_sum + YearProposed, family = "poisson", data = sc_subs_yearly)
summary(model14)
logLik(model14)
AIC(model14)

model15 <- glm(low_strength_proposed ~ us_rus_dist + meeting_number + veto_sum + YearProposed, family = "poisson", data = sc_subs_yearly)
summary(model15)
logLik(model15)
AIC(model15)

##--Table 13 appendix: SBs Created, Ordered logit

#NOTE: Need to download meeting level data for these regressions
meeting_level <- read.csv(file = "unsc_sb-meetinglevel_final.csv")

#estimate ordered models 
model1.ordered <- polr(as.factor(Strength) ~ p5_avg + meeting_n + veto_sum, data = meeting_level, Hess=TRUE)
summary(model1.ordered)

model2.ordered <- polr(as.factor(Strength) ~ p5_avg + meeting_n + veto_sum + log(pk_forces) + cold_war, data = meeting_level, Hess=TRUE)
summary(model2.ordered)

model3.ordered <- polr(as.factor(Strength) ~ IdealPointDistance + meeting_n + veto_sum, data = meeting_level, Hess=TRUE)
summary(model3.ordered)

model4.ordered <- polr(as.factor(Strength) ~ IdealPointDistance + meeting_n + veto_sum + log(pk_forces) + cold_war, data = meeting_level, Hess=TRUE)
summary(model4.ordered)

##get p-values for ordered models
ctable1 <- coef(summary(model1.ordered))
p1 <- pnorm(abs(ctable1[,"t value"]), lower.tail = FALSE*2)
ctable1 <- cbind(ctable1, "p value" = p1)
ctable1

ctable2 <- coef(summary(model2.ordered))
p2 <- pnorm(abs(ctable2[,"t value"]), lower.tail = FALSE*2)
ctable2 <- cbind(ctable2, "p value" = p2)
ctable2

ctable3 <- coef(summary(model3.ordered))
p3 <- pnorm(abs(ctable3[,"t value"]), lower.tail = FALSE*2)
ctable3 <- cbind(ctable3, "p value" = p3)
ctable3

ctable4 <- coef(summary(model4.ordered))
p4 <- pnorm(abs(ctable4[,"t value"]), lower.tail = FALSE*2)
ctable4 <- cbind(ctable4, "p value" = p4)
ctable4

